#This file produces the responsiveness measure used in section 4 of the Cambridge Element
#Run on MacOS Monterey 12.6
#R version 4.4.0 (2024-04-24)

#Load packages
rm(list = ls()) #clear environment
library(tidyverse)
library(coda)
library(car)
library(dplyr)
library(ggplot2)
library(rjags)


######################################
#Import and organise divisions data
######################################

#Import MP names 
mpnames <- read.delim("./Other data sources/mps-2017.txt", header = TRUE, sep = "\t")
mpnames <- unite(mpnames, mpname, c(firstname, surname), remove=TRUE, sep = " ") #merge first name and surname
mpnames <- mpnames[,c(1,2)] #keep only MP name and id
mpnames$mpname <- gsub("O&#39;", "O'", mpnames$mpname) #fix encoding issue 

#Import roll call matrix
#transform into data frame with votenumbers as row names, mp names as column names 
rollmatrix <- read.delim("./Other data sources/votematrix-2017.txt", header = TRUE, sep = "\t")
rollmatrix$Bill <- gsub("&#8217;", "'", rollmatrix$Bill) #fix encoding issue 
rollmatrix$date <- as.character(rollmatrix$date)
rollmatrix$voteno <- as.character(rollmatrix$voteno)

#Select Brexit-related bills
#see document called BrexitBillSelect.docx to understand which bills are deemed relevant.
inxvoteno <- grep("293|307|308|310|311|312|331|332|333|345|346|347|354|357|358|359|360|361|362|363|364|374|386|387|388|389|390|391|392|393|394|395|397|398|399|400|404|405|406|407|408|409|439|440|441|442|\\<3\\>|\\<4\\>|\\<5\\>|\\<11\\>", rollmatrix$voteno, fixed = F)
rollmatrix <- rollmatrix[inxvoteno,]
rollmatrix <- rollmatrix[rollmatrix$date > as.Date("2019-01-01"),]
rownames(rollmatrix) <- rollmatrix$voteno # realign row names
rolls <- rollmatrix[,-c(1,2,3,4)] #drop some columns (rowid, date, voteno, bill)
rolls <- subset(rolls, select=-X) #remove this empty column named 'X'

#Transform vote matrix
colnames(rolls) <- gsub("mpid", "", colnames(rolls)) #remove "mpid" from variables
rolls <- t(rolls) #transpose matrix

#Replace voting codes
rolls <- recode(rolls, "2=10") #yes
rolls <- recode(rolls, "1=10") #yes
rolls <- recode(rolls, "-9=0") #NA
rolls <- recode(rolls, "3=0") #NA
rolls <- recode(rolls, "4=20") #no
rolls <- recode(rolls, "5=20") #no

#Join MPs' names
rolls <- cbind(rolls, mpid = unlist(rownames(rolls)))
rolls <- merge(mpnames, rolls, by = "mpid", all.y = TRUE) #match mp names with mp ids 
rolls <- cbind(rolls[,1:2], mutate_all(rolls[,3:52], function(x) if(is.character(x)) as.numeric(as.character(x)) else x)) #change votes back to numeric

#Collapse data by MP (for MPs who switched party affiliation)
rolls <- rolls %>%
  group_by(mpname) %>%
  summarise_all(list(max)) #summarize voting records by mp (deprecated function warning but works just fine)

#Recode again to a 0-1-NA format necessary for IRT models 
rolls <- replace(rolls, rolls == 0, NA) #NA
rolls <- replace(rolls, rolls == 10, 1) #YES
rolls <- replace(rolls, rolls == 20, 0) #NO

#Get rid of MPs with huge gaps in voting record
rolls$na_count_p1 <- apply(rolls[,11:52], 1, function(x) sum(is.na(x))) #no of missing votes in period 1
rolls$na_count_p2 <- apply(rolls[,3:10], 1, function(x) sum(is.na(x))) #no of missing votes in period 2
rolls <- rolls[rolls$na_count_p1 <= 21 & rolls$na_count_p2 <= 4,] #exclude if more than half of votes in either period missing

#Identify anchors
rolls$anchorMP <- NA
rolls$anchorMP[rolls$mpname == "Steven Baker"] <- 1
rolls$anchorMP[rolls$mpname == "David Lammy"] <- 2
rolls <- rolls[order(rolls$anchorMP),-c(53)] #put anchor MPs to the top

#Restrict to Ys
rolls.jags <- rolls[,3:52]
rolls.jags <- as.matrix(rolls.jags)


###########################################
#Dynamic IRT Model with Quadratic Effects
###########################################

#JAGS code for model with fixed evolution parameter
dyn_irt_model_quadratic <- 'model{
#loop through actors
for(i in 1:nactors){
#loop through votes
for(j in 1:nvotes){
Y[i,j] ~ dbern(mu[i,j])
logit(mu[i,j]) <- alpha[j] + beta[j] * theta[actor[i], period[j]] + beta2[j] * theta[actor[i], period[j]] * theta[actor[i], period[j]]
}
change[i] <- theta[i,2] - theta[i,1]
}
#set normal priors on betas and alphas
for(j in 1:nvotes){
beta[j] ~ dnorm(0, 1)
beta2[j] ~ dnorm(0, 1)
alpha[j] ~ dnorm(0, 1)
}
#set priors on thetas (fix space with two actors)
theta[1, 1] ~ dnorm(1, 1)
theta[2, 1] ~ dnorm(-1, 1)
for(t in 2:nperiods){
theta[1, t] ~ dnorm(theta[1, t-1], tau.evol)
theta[2, t] ~ dnorm(theta[2, t-1], tau.evol)
}
for(c in 3:nactors){
theta[c, 1] ~ dnorm(0, 1)
for(t in 2:nperiods){
theta[c, t] ~ dnorm(theta[c, t-1], tau.evol)
}
}
#set priors on evolution parameters for thetas
tau.evol <- 1
}
'

#Set parameters
jags.data.irt <- list(
  Y = rolls.jags,
  nactors = 630,
  nvotes = 50,
  nperiods = 2,
  period = c(2, 2, 2, 2, 2, 2, 2, 2, rep(1, 42)),
  actor = as.numeric(c(1:630))
)

#Fit JAGS model
model.jags.irt <- jags.model(file=textConnection(dyn_irt_model_quadratic), data = jags.data.irt, n.chains = 1, inits=list(.RNG.name="base::Wichmann-Hill",.RNG.seed=2355328))
posterior.coda.irt <- coda.samples(model.jags.irt, c("theta","beta","beta2","alpha", "tau.evol", "change"), n.iter=10000, thin = 10, progress.bar="text")

#Extract debate effects (betas)
sum.jags.irt=summary(posterior.coda.irt)[[1]][,1]
sum.jags.irt.beta=data.frame(means=sum.jags.irt[grepl(names(sum.jags.irt),pattern="beta")])
sum.jags.irt.beta$debate=as.numeric(gsub("\\D","",row.names(sum.jags.irt.beta)))

#Extract actor positions (thetas)
sum.jags.irt.theta=data.frame(means_period1_sq=sum.jags.irt[grepl(names(sum.jags.irt),pattern="theta")&
                                                           grepl(names(sum.jags.irt),pattern=",1")],
                              means_period2_sq=sum.jags.irt[grepl(names(sum.jags.irt),pattern="theta")&
                                                           grepl(names(sum.jags.irt),pattern=",2")],
                              change_sq=sum.jags.irt[grepl(names(sum.jags.irt),pattern="change")])

#Get 95% credible intervals
brexit.irt=HPDinterval(posterior.coda.irt[[1]],prob=.95)

#Save CIs in data frame
brexit.irt.thetas=data.frame(lower_period1_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "theta")&grepl(row.names(brexit.irt),pattern = ",1"),"lower"],
                          upper_period1_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "theta")&grepl(row.names(brexit.irt),pattern = ",1"),"upper"],
                          lower_period2_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "theta")&grepl(row.names(brexit.irt),pattern = ",2"),"lower"],
                          upper_period2_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "theta")&grepl(row.names(brexit.irt),pattern = ",2"),"upper"],
                          lower_change_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "change"),"lower"],
                          upper_change_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "change"),"upper"])

#Add MP names and actor numbers
brexit.irt.thetas$MP=rolls$mpname
brexit.irt.thetas$MPnumber=1:nrow(rolls)
sum.jags.irt.theta$MP=rolls$mpname
sum.jags.irt.theta$MPnumber=1:nrow(rolls)

#Join CIs and means
brexit.irt.thetas_sq=left_join(brexit.irt.thetas, sum.jags.irt.theta)

#Assess convergence of model with Geweke diagnostic
geweke_irt <- unlist(geweke.diag(posterior.coda.irt))
geweke_irt.dat=data.frame(gew=geweke_irt[!grepl("change",names(geweke_irt))&!grepl("frac",names(geweke_irt))],
                          names=names(geweke_irt[!grepl("change",names(geweke_irt))&!grepl("frac",names(geweke_irt))]))
geweke_irt.dat=geweke_irt.dat%>%filter(!is.na(gew))

#Assess how many parameters fall in 2/-2 range
prop.table(table(geweke_irt.dat$gew>2|geweke_irt.dat$gew< -2))
table(geweke_irt.dat$gew>2|geweke_irt.dat$gew< -2) #looks as expected - good convergence

#Plot Geweke diagnostics
geweke_irt_graph=ggplot(geweke_irt.dat) + 
  geom_histogram(binwidth=0.2,aes(x=gew,y=..density..), position="identity") + 
  geom_density(aes(x=gew,y=..density..))+
  labs(x="Geweke Statistics",y="Density")+
  theme_bw()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        axis.text=element_text(size=12,colour="black",family="Verdana"),
        axis.title=element_text(size=12,family="Verdana"),
        legend.position="none")+
  scale_x_continuous(limits = c(-4,4))
geweke_irt_graph #looks like convergence

#Rename data
MP_data <- brexit.irt.thetas_sq

#Export data
write.csv(MP_data, "./Created auxiliary datasets/irt_models_MPs.csv")


